% [epoch badChan]=bciGetEpoch(bci,dat)
% creates epochs of data segment according to parameter setting in
% bciGetParam (which is recommended to be located in the subject's folder)

function [epoch badChan]=bciGetP3Epochs(bci,dat,prm,onsets,N)

doDebug = true;

% cast double precision, otherwise absurd results with filtfilt and maybe
% others
dat = double(dat);

% do baseline correction
if bci.init.nBaseline>=1
    baseline = repmat(mean(dat(:,1:bci.init.nBaseline),2),1,bci.init.nEpochLength);
    dat = dat - baseline;
end


% detrend time series
if bci.param.detrend,
    dat=detrend(dat')';
end

% find artifacts
badChan = sum(abs(dat(:,bci.init.intervalOfInterest))>bci.param.artifactThresh,2) > 0;

% remove common average reference
if ~isempty(bci.refChanCAR)
    ecog.data = dat;
    ecog.selectedChannels = bci.refChanCAR;%1:size(dat,1); % exclude bad channels of this single epoch???
    ecog=ecogRemoveCommonAverageReference(ecog,'equal');
%     ecog=ecogRemoveCommonAverageReference(ecog,'projected');
    dat = ecog.data;
end

% spatial filter
if bci.param.spatialFilter>0,
    dat = (dat'*bci.init.spatialFilter)';
end

% bandpass filter, important: cast to double
if length(bci.param.bandFreqs)>0,
    dat = filtfilt(bci.init.filtCoeffB,bci.init.filtCoeffA,dat')';
end

if  bci.param.specMethod==0,
    % timeseries data
    if ~isempty(bci.param.resampFreq),
        % resample
        dat=resample(dat', floor(bci.param.resampFreq), floor(bci.srate))';
        intOfI=bci.init.resampIntOfInterest;
    else
        intOfI=bci.init.intervalOfInterest;
    end
    dat=dat(:,intOfI);
else
    error('bci.param.specMethod not defined');
end

% select the flash specific epochs
% initialize
flashWinLength = round(bci.param.resampFreq*bci.paradigm.flashWin); % samples per epoch
idxTmpl = 0:flashWinLength-1;

% 
flashSequenceMat = prm.flashSequenceMat(:,prm.tobjIdx,N);
flashSequenceMat(sum(flashSequenceMat,2)==0,:)=[];

rtime = linspace(0,size(dat,2)/bci.param.resampFreq,size(dat,2));
realFlashOnsetTimes = onsets/bci.srate;
isValidPattern = [ones(1,size(dat,2)),zeros(1,flashWinLength)];

if bci.paradigm.doAverageFlashEvents,        
    targetEp = zeros(size(dat,1),flashWinLength,length(prm.tobjIdx));
    nonTargetEp = zeros(size(dat,1),flashWinLength,length(prm.tobjIdx));
    targetEpochs = zeros(size(dat,1),flashWinLength,sum(flashSequenceMat(:,1)));
    nonTargetEpochs = zeros(size(dat,1),flashWinLength,sum(~flashSequenceMat(:,1)));
    % ensure involvement of all target flashes in average
    nTargetInMean = zeros(1,flashWinLength);
    % get the mean for each target/nontarget
    for target = 1:length(prm.tobjIdx),
        targetCount = 0;
        nonTargetCount = 0;
        for tp = 1:length(realFlashOnsetTimes);
            onset = find(rtime>realFlashOnsetTimes(tp),1);
            if isempty(onset), continue; end
            isValid = isValidPattern(onset+idxTmpl)>0;
            % cut epochs
            if flashSequenceMat(tp,target),
                targetCount = targetCount+1;
                targetEpochs(:,isValid,targetCount) = dat(:,onset+idxTmpl(isValid));
                nTargetInMean = nTargetInMean + double(isValid); % ignore outlying samples in average
            elseif all(isValid==1) % include only fully available non-target sequences
                nonTargetCount = nonTargetCount + 1;
                nonTargetEpochs(:,:,nonTargetCount) = dat(:,onset+idxTmpl);
            end
        end
        % get the mean
        targetEp(:,:,target) = sum(targetEpochs,3)./repmat(nTargetInMean,size(dat,1),1);
        nonTargetEp(:,:,target) = mean(nonTargetEpochs,3);
    end
    % reshape the epochs to the form:
    % Target1features,Target2features,..., NonTarget1features,NonTarget2features,...
    if isfield(prm,'positionMat'),
        reOrderIdx = prm.positionMat(:,N);
    else
        reOrderIdx = 1:prm.nObjects;
    end
    epoch = cat(2,reshape(targetEp(:,:,reOrderIdx),[size(targetEp,1),size(targetEp,2)*size(targetEp,3)]),...
        reshape(nonTargetEp(:,:,reOrderIdx),[size(nonTargetEp,1),size(nonTargetEp,2)*size(nonTargetEp,3)]));
else
    % sort target flash occurence by object ID
    if any(sum(flashSequenceMat,2)~=1),
        warning('Non-uniqe number of target flashes. Single flash event classification not supported.');
    end
    flashesPerTarget = unique(sum(flashSequenceMat,1));
    [flashObj,dummy]=find(flashSequenceMat'==1);
    sortIdx = zeros(1,flashesPerTarget*length(prm.tobjIdx));
    for target = 1: length(prm.tobjIdx),
        sortIdx((target-1)*flashesPerTarget+1:target*flashesPerTarget)=find(flashObj==target);
    end
    epoch = zeros(size(dat,1),flashWinLength*size(flashSequenceMat,1));
    curEpoch = zeros(size(dat,1),flashWinLength);
    lastValid = zeros(1,length(prm.tobjIdx));
    for tp = 1:length(sortIdx);
        onset = find(rtime>realFlashOnsetTimes(sortIdx(tp)),1);
        if isempty(onset), 
            warning('Empty P3 epoch detected.');
            isValid = false(1,length(idxTmpl));         
        else
            isValid = isValidPattern(onset+idxTmpl)>0;
        end
        curEpoch(:,isValid)=dat(:,onset+idxTmpl(isValid));
        if any(~isValid), % mix a former epoch with the current one
            tp2=lastValid(flashObj(sortIdx(tp)));
            epIdx = (tp2-1)*flashWinLength+1:tp2*flashWinLength;
            mixingWeight = ones(size(epoch,1),1)*[linspace(1,0,sum(isValid))*find(sum(isValid)),zeros(1,sum(~isValid))];
            curEpoch=curEpoch.*mixingWeight + epoch(:,epIdx).*(1-mixingWeight);
        else
            lastValid(flashObj(sortIdx(tp)))=sortIdx(tp);
        end
        epoch(:,(tp-1)*flashWinLength+1:tp*flashWinLength)=curEpoch;
    end    
end
